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We present solutions of coupled particle-field evolution in classical U{1) and SU{2) gauge theories 
in real time on three-dimensional lattices. Our simulations are performed in a regime of extreme 
C^ , anisotropy of the momentum distribution of hard particles where back-reaction is important. We 

find qualitatively different behavior for the two theories when the field strength is high enough that 
non-Abelian self-interactions matter for SU{2). It appears that the energy drained by a Weibel-like 
plasma instability from the particles does not build up exponentially in soft transverse magnetic 
fields but instead returns, isotropically, to the hard scale via a rapid avalanche into the ultraviolet. 
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I. INTRODUCTION 

One of the most important open questions emerging from the Relativistic Heavy-Ion 
CoUider (RHIC) program is how quickly, and by what process, the stress-energy tensor of 
the high-density QCD matter produced in the central rapidity region becomes isotropic. 
This is important, for example, to determine the maximal temperature achieved in such 
collisions, which will also be pursued in the near future at even higher energies at the CERN 
Large Hadron Collider (LHC). One of the chief obstacles to thermalization in ultrarelativistic 
heavy-ion collisions is the rapid longitudinal expansion of the central rapidity region. If the 
matter expands too quickly then there will not be sufficient time for its constituents to 
interact and thermalize. In the absence of interactions, the longitudinal expansion causes 
the system to become much colder in the longitudinal direction than in the transverse 
directions, corresponding to (p|) <^ {p^) in the local rest frame. One can then ask how long 
it would take for interactions to restore isotropy in momentum space. 

High-energy heavy-ion collisions release a large amount of partons from the wavefunctions 
of the colliding nuclei. Partons with very large transverse momenta originate from high-Q^ 
hard interactions, while partons with transverse momenta below the so-called saturation 
momentum Qs (given by the square root of the color charge density per unit area in the 
incomingnuclei) are much more abundant and are better viewed as a classical non-Abelian 
field pj, y, l3|. At high energies, Qs 3> Aqcd sets a semi-hard scale which allows for a 
weak-coupling treatment of the early-stage dynamics. However, due to the large occupation 
number of the soft modes below Qg (the classical color field has strength F^^ ~ ^1 9)-i the 
problem is nevertheless non-perturbative. 

The evolution following the initial impact was among the questions which the so-called 
"bottom-up scenario" jj] attempted to answer. For the first time, it addressed the dynamics 
of soft modes ("fields") with momenta much below Qs coupled to the hard modes ("par- 
ticles") with momenta on the order of Qs and above. However, it has emerged recently 
that one of the assumptions made in this model is not correct. The debate centers around 
the very first stage of the bottom-up scenario in which it was assumed that (a) collisions 
between the high-momentum (or hard) modes were the driving force behind isotropization 
and that (b) the low-momentum (or soft) fields act only to screen the electric interaction. In 
doing so, the bottom-up scenario implicitly assumed that the underlying soft gauge modes 
behaved the same in an anisotropic plasma as in an isotropic one. However, it turns out 
that in anisotropic plasmas, within the so-called "hard loop approximation" (see below), the 
most important collective mode corresponds to an instability to transverse magnetic field 
fluctuations [a, la, l7|. Recent works have shown that the presence of these instabilities is 
generic for distributions which possess a momentum-space anisotropy [sl, [9, [lO| and have 



obtained the full hard-loop (HL) action in the presence of an anisotropy [ll|]. Another im- 
portant development has been the demonstration that such instabilities exist in numerical 



solutions to classical Yang- Mills fields in an expanding geometry [12, ll3|, [14 



In the last year there have been significant advances in the understanding of non-Abelian 



soft-field dynamics in anisotropic plasmas within the HL framework |15l . 1 161 . |17I |. The HL 
framework is equivalent to the collisionless Vlasov theory of eikonalized hard particles, i.e. 
the particle trajectories are assumed to be unaffected by the induced background field. It 
is strictly applicable only when there is a very large scale separation between the soft and 



hard momentum scales.^ Even with these simphfying assumptions, HL dynamics for self- 
interacting gauge fields is non-trivial due to the presence of non-linear interactions which 
can act to regulate unstable growth. These non-linear interactions become important when 
the vector potential amplitude is on the order of Anon-Abeiian ~ Vsl 9 ~ ^ffhPh, where ph is 
the characteristic momentum of the hard particles, fh is the angle-averaged hard occupancy, 
and ps ~ g\f7hPh is the characteristic soft momentum of the fields.^ In QED there is no 
such complication and the fields grow exponentially until AAbeUan ~ Ph/ 9 at which point the 
hard particles undergo large-angle deflections by the soft background field invalidating the 
assumptions underpinning the hard-loop effective action. 

Recent numerical studies of HL gauge dynamics for SU(2) gauge theory indicate that for 
moderate anisotropics the gauge field dynamics changes from exponential field growth indica- 
tive of a conventional Abelian plasma instability to linear growth when the vector potential 
amplitude reaches the non- Abelian scale, Anon-Abciian ~ Ph (l5|, ll6|. This linear growth 
regime is characterized by a cascade of the energy pumped into the soft modes by the insta- 



bility to higher-momentum plasmon-like modes 28|, |29|. These results indicate that there 
is a fundamental difference between Abelian and non-Abelian plasma instabilites. However, 
even with this new understanding the HL framework relies on the existence of a very large 
separation between the hard and soft momentum scales by design. One would like to know 
what happens when the back-reaction on the hard modes is not completely negligible, which 
is probably the situation faced in real experiments at finite energies. In this case one is nat- 
urally led to consider instead the Wong- Yang- Mills (WYM) equations. These equations can 



be shown to reproduce the HL effective action in the weak-field approximation [30|, l31|, l32 
however, when solved fully they go beyond the HL approximation. 

In this paper we present first real-time three-dimensional lattice results obtained by so- 
lution of the WYM equations for coupled U(l) and SU(2) particle-field systems. For SU(2) 
gauge group, the WYM equations describe the propagation of classical colored particles in a 
colored background field in a framework that self-consistently includes the deflection of the 
hard particles by the soft fields. In practice, the equations are solved using the test-particle 
method and we present an improved numerical algorithm which uses current smearing to 
enable us to obtain convergent results with much fewer test particles than the straightfor- 
ward "nearest grid point" (point particle) method. We focus here on simulations with large 
anisotropics of the particles and, for SU(2), on initial field strengths beyond the non-Abelian 
scale. 

Our numerical results agree in some qualitative aspects with the HL calculations in that 
we see a transfer of energy from particles to the soft fields which saturates at late times. 
However, for the simulations with large anisotropy presented here, the saturation seems 
to stem from an "avalanche" of energy which has been dumped in soft modes to higher 
momentum modes. Due to finite lattice spacing this avalanche stops at the highest momen- 
tum modes of our lattice at which time the rapid growth of the field energy density ends. 
As the lattice spacing in our simulations is decreased we find that the amplitude at which 
the field energy density saturates increases, indicating that the mechanism for saturation 
is the population of hard lattice modes. This is the main difference to the cascade which 



^ Of course, in the hard-loop approach very small angle deflections, 9 ^ g, are included self-consistently. 
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We will not attempt here to use the s pec ific fields obtained from the Color Glass Condensate wave 
functions [l8|, UM, l20|, l2l|, [22, l23|, l2J, |25|, l26|, l27| as initial conditions. Hence, we label the hard particle 



momentum scale more generically as ph rather than Qs 



occurs in hard-loop simulations with moderate anisotropics [29|, where field modes near 
the hard particle scale do not get populated until parametrically large times. The behav- 
ior seen here might be related to the chaoticity of three-dimensional classical Yang-Mills 



fields (33|, |3J, |35| . Additional studies are necessary, however, before firm conclusions in this 
regard can be drawn. 

In Sec. HI] we present the Wong- Yang-Mills equations and the basics of the test particle 
method. In Sec. IIIII we discuss the numerical methods employed. In Sec. HV] we present 
results for U(l) and SU(2) gauge theories. In Sec. |V] we discuss the implications of our 
results and list open questions which remain. 



II. WONG- YANG-MILLS EQUATIONS 



In this paper we present solutions of the classical Vlasov transport equation for hard 
gluons with non-Abelian color charge g" in the collisionless approximation 36l. l37l|. 



p''[d, - gq^F^B; - gU,AU'^d,a]f{x,p, q) = 



(2.1) 



Here, /(t, x,p, q"-) denotes the single-particle phase space distribution function and g is the 
gauge coupling, which can take arbitrary values at the classical level. 

The Vlasov equation is coupled self-consistently to the Yang-Mills equation for the soft 
gluon fields, 

9 [^dqqv''f{t,x,p,q), (2.2) 
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with v^ = (1, v/p). These equations reproduce the "hard thermal loop" effective action 
near equilibrium J30|, ImI, l32]- However, the full classical transport theory (12. 1112. 2p also 
reproduces some higher n-point vertices of the dimensionally reduced effective action for 
static gluons |38| beyond the hard-loop approximation. The back-reaction of the long- 
wavelength fields on the hard particles ("bending" of their trajectories) is, of course, taken 
into account. This is important for understanding particle dynamics in strong fields or for 
extremely anisotropic particle momentum distributions. 

Eq. (12.11) can be solved numerically by replacing the continuous single-particle distribution 
f{x,p, q) by a large number of test particles: 

/(^, P^<i) = ^Y. ^(=^ - ^'(^)) (2^)'^(P - P^ W) ^(^" - ^"(^)) ' (2-3) 

-/Vtpst 



where Xi(t), Pi{t) and qfit) are the coordinates of an individual test particle and Attest 
denotes the number of test-particles per physical particle. The Ansatz (12. 3p leads to Wong's 
equations |36|, |37| 



dxi 
dt 


= Vi, 


dpi 
df 


= gqtiE'^ + v.xB'^), 


dqi 
dt 


= tg v'^ [A^, Qi], 


jau 


= ThT^Qtv^Six-x^it)) 



(2.4) 

(2.5) 

(2.6) 
(2.7) 



'tost 



for the i-th test particle. The time evolution of the Yang-Mills field can be followed by the 
standard Hamiltonian method |39|] in A^ = gauge. Specific numerical algorithms to solve 
the classical field equations coupled to particles will be presented in the next section. 
We use the following dimensionless lattice variables: 



9(^ jpa ^ _ a^ ^a _ 1 a .r _^t 



El = '—E% Pl = ^P, Ql = ,^q\ iVtcst,L = ^. (2.8) 

Our lattice Hamiltonian^ is then 

^^ = \Y.^Vl + \Y. (^-^ - ^^ Tr Uu) + -^ Y. IP^-.^I' (2.9) 

which is related to the physical Hamiltonian hy H = {4/g'^a) Hi. To convert lattice vari- 
ables to physical units we fix the length of the lattice to L = 5 fm which then determines 
the physical scale for the lattice spacing a. All other dimensionful quantities can then be 
converted to physical units using eqs. (12. 8p . In other words, simulations on lattices with 
increasing number of sites correspond to a fixed infrared cutoff, while the ultraviolet cutoff 
increases in proportion to the number of sites. In lattice units, ph,L should not be much 
smaller than 1; otherwise, the particle and field modes would overlap significantly. 

For most of the results presented here the initial phase-space distribution of hard gluons 
is taken to be a (parity-invariant) planar momentum distribution: 

f{p) = ng i^—j 5{p,) eM-PT/Ph) , (2.10) 



with pt = v/Pi + pI aiid Hg the number density of hard gluons (summed over two helicities 
and, for SU{Nc) gauge group, also over A''^ — 1 color states). This represents a quasi-thermal 
distribution in two dimensions, with average momentum equal to ph- At the initial time 
t = we randomly sample Np = Attest Ug a^ particles from the distribution (12.101) at each cell 
of the lattice.^ 

We shall also present data for distributions with more moderate or vanishing anisotropy 
for numerical checks, where indicated. Instead of the number density Ug of hard gluons it is 
also customary to quote the assumed value for the mass scale 

ml-9^NjjX^~9^NA. (2.11) 



(27r)3 \p\ pf^ 

(For the particular distribution (I2.10p the latter is actually an equality.) This quantity sets 
the scale for the growth rate of unstable field modes in the linear approximation. 

The initial field amplitudes are sampled from a Gaussian distribution with a width tuned 
to a given initial energy density: 



{A'^{x)A%y)) = -^5,,5-'5{x - y) , (2.12) 



^ We employ the compact lattice formulation even for C/(l) gauge group. 

^ When Np is not very large it is useful to ensure explicitly that the sum of particle momenta in each cell 
vanishes, for example by adjusting the momentum of the last particle accordingly. 



and E = —A = 0. Gauss's law then implies that the local charge density at time t = 
vanishes. We ensure that any particular initial condition satisfies exact local charge neu- 
trality. In our U{1) simulations we monitor the time evolution of Gauss's law as a check 
of the numerical accuracy. Our charge smearing algorithm for SU{2), on the other hand, 
explicitly exploits (covariant) current conservation and hence Gauss's law is satisfied exactly 
by construction. For f/(l) simulations the results of single simulation runs are shown but 
for our SU{2) simulations observables are averaged over several initial field and particle 
configurations. Where shown error bars indicate the root-mean-square variation over this 
ensemble. 



III. NUMERICAL METHODS 

A. NGP method for non-Abelian gauge theories 

We first summarize the method developed in Ref. J40|, |4l| where zero-order weighting is 
used (i.e., point particles). In that approach, one simply counts the number of particles 
within a cell to determine the charge density, which is why it is called the nearest-grid-point 
(NGP) method. A current is generated only when a particle crosses a cell boundary from i 
to i + 1. This induces a current on that link which is J{i) = Q6(t — tcross)/^tcst- In order to 
satisfy the lattice covariant continuity equation, 

A = X] Ulii - x)J^{i - x)U^{i - x) - J^{i) , (3.1) 

X 

the charge must be parallel transported: 

Q{i + l) = Ul{i)Q{i)U,{i) . (3.2) 

The rotation of the particle's momentum by the magnetic field can be updated continuously 
but its magnitude changes only when it crosses a link. The momentum "kick" can be 
determined by energy conservation: 

l_ I . -i'test j-,2 I I , -''test/,-, tn2 /o o\ 

IPinil + ^— -^ini = IPfinI + ^"(^ini " J) ■ (3.3) 

Substituting J = Q/Nt^st and assuming that the particle crosses the cell boundary in x- 
direction, one obtains 

IPfinI = IPinil - ^x.miQ + QV{2Ntest) • (3.4) 

This way, the total energy 

lattice i 

is conserved. 

In prin ciple, the NGP method also requires time ordering of the link crossings by par- 



ticles J40|, |41i in every time step, which becomes quite time consuming when the number 
of test-particles is large. However, we have found that at least for our present purposes 
time-ordering does not affect the results significantly and so will be ignored in what follows. 
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This method was apphed in 4^, |43|, |4^ to study a simphfied situation with fields and par- 
ticles living on a one-dimensional spatial lattice. Physically, this corresponds to a transver- 
sally homogeneous system (translational invariance in x and y directions). For such ld-3v 
simulations ^ the current can be made sufficiently smooth by employing enough test-particles. 
Ref. |42|, |43|, |4J] found that the currents were sufficiently smooth when the number of par- 



ticles per lattice site, Np 



Ntest rig a-^ 



was on the order of a few hundred to a thousand. 



50 

^ 45 

CM 

I 35 
I 30 
I 25 
u> 20 
c 15 

LU 

10 



^-B2/2 

8^/2 

■E2/2 

£2/2 



U(1)16^Np=1000NGP 
Isotropic 




10 



20 



30 



50 

^ 45 

CM 

,-.40 
I 35 
I 30 
I 25 

S; 20 

c 15 

LU 



10 







^-6^/2 


U(1) 16^ Np=50 Smearing 


: ° B^/2 


Isotropic 


: ■-£2/2 




^ - £2/2 


, . .-lo 


:. ,Jl4-,Q^SAff»- 


AS/S£«^/ y^ y ^PxT* -"/^ '/^^.^^'''sJfS^jf 


s^^yCJ^" 


^St'-^^x^. a Ar*\/* * V/^^'^'* *■" ■ 




" ' 


- 





10 



20 



30 



mj 



mj 



FIG. 1: Time evolution of the magnetic and electric energy densities for U{\) gauge group and 
isotropic particle momentum distribution using the nearest-grid-point method (left) and smeared- 
particles (right). 

Full 3d-3v simulations, however, require improved algorithms since 3d lattices have many 
more sites and the NGP method becomes computationally too expensive. In Fig. [1] (left), 
the time evolution of the field energy densities on a 16'^ lattice with Np = 1000 test-particles 
per site are shown. Despite the large number of particles, there seems to be an anomalous 
damping of the fields with the NGP method. We naively expect that to achieve the same 
accuracy as for a one-dimensional lattice, Np would have to increase exponentially with 
the number of spatial dimensions: A'^'' ~ {Np'^)^. This makes multi-dimensional NGP 
simulations practically impossible and forces one to consider methods for smoothing particle 
densities and currents. 



B. Current smearing in electromagnetic PIC simulations 



In Abelian plasmas it is common to employ smoothed currents for particle-in-cell (PIC) 
simulations |45|, |46| in order to avoid numerical noise. The charge density at site {i,j, k) is 
obtained by smearing ^ 



Nr, 



p{i,j,k) = e y^ Si{Xn)Sj{yn)Sk{Zn) 



(3.6) 



n=l 



^ Fields and charge densities fluctuate only in z-direction but particle velocities can point in any direction. 
^ We assume Ax = Ay = Az = 1. 



Here, S* is a form factor. For example, the first-order form factor (shape-factor) is 

qUt\ ^ / 1 - 1^ - ^1 for 1^ - z| < 1, 
■'^^^ \0 for|^-z|>l 



(3.7) 



However, it is well known that if we naively define the current density at each grid point as 



Ja{hj,k) = e^Va,nSi{Xn)Sj{yn)Sk{Zn), a = x,y, 



(31 



n=\ 



then the continuity equation is not satisfied exactly [45|, |46|, requiring us to correct the 
electric field by solving the Poisson equation. 

To avoid solving Poisson's equation in each time step, several numerical techniques for 
solving the continuity equation have been developed [47|, |48|, |49|, l50|, ImI]- For simplicity, we 
consider the 2d problem to explain how to satisfy local charge conservation in electromag- 
netism. We consider a particle moving from (xi,yi) to {x2,y2) during a time step At, i.e. 
2^2 = a;i -|- VxAt, 1/2 = 2/1 + VyAt. We also restrict ourselves to the case in which (0:2,2/2) 
belongs to the same lattice site {i,j)- Charge conservation of the assigned current densities 
is realized with the procedure described by Eastwood 47l |. 



Jyihj) 



Fy{l-W,), 



F.Wy, 



Jy{l+l,j)=FyW, 



y" X 1 



where {Fx,Fy) = F represents the charge flux 

_ X2-X1 
^""^ At ' 



rj. 2/2-2/1 

^V = e 7 

^ At 



(3.9) 



(3.10) 



W represents the first-order shape-factor corresponding to the linear weighting function 
defined at the midpoint between the starting point (xi,2/i) and the end point (3:2,2/2) 



W., 



Xi + X2 



I , 



W„ 



2/1 + 2/2 



2 ' ' 2 

One can easily show that the lattice continuity equation 



Jxihj) - Jx{i - 1, j) + Jy{id) - Jyihj - 1) 



P[hJ) 



J ■ 



p*+^*(z,j; 



At 



(3.11) 



(3.12) 



is satisfied with the definition (13. 9p of the current. When (0:2,2/2) is a different lattice site 
than the original («, j) one, we use the 'Zigzag scheme' developed in Ref. 5l|. Much of the 
numerical noise is eliminated by such linear interpolation. Note that our current may in 
general be expressed as 



N„ 



Jaihj^k) = e^Va,nS^{xn)SJ{yn)Sl{zn), a = x,y,z, aj^j,k, (3.13) 



n=l 



which differs from (13.81) since 5*° = dS\/dXa] 5*° is the "form factor" of the NGP method, 
i.e. 5"° = 1 if the particle is in the corresponding lattice site, otherwise it is zero. 
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Finally, we note that the electromagnetic forces should be smeared in a similar way when 



a particle momentum is updated [47|, |48|]. Consider the time derivative of the total energy: 



^ = _^^.j+ J2 q^E{x,{t))-v,. (3.14) 

lattice particles 

It is clear that the interpolation function for E should be the same as that for J in order to 
achieve good energy conservation in the simulation. The electric field E[x) at the particle 
position X is then obtained from 

E^{x) = ^ SlSlS}^E^{z,j,k) , (3.15) 

lattice 

while the magnetic field is given by 

B^{x) = J2 SlSlS?^BS,m . (3.16) 

lattice 

This is motivated by the relations E = —A and B = 'V x A. The momentum update itself 
is explained in appendix \M 

In Fig. [T] (right) we show the U{1) energy densities resulting from an isotropic run using 
Np = 50 smeared test-particles per site. As can be seen from this figure the smeared-particle 
result does not suffer from the anomalous damping seen with the NGP method. More 
importantly, this result was obtained using many fewer particles than would be required to 
obtain a stable result with the NGP method. Therefore the smeared-particle method makes 
three-dimensional particle-in-cell (PIC) simulations possible in practice. 

Fig. [2] shows the f/(l) field energy densities as a function of time for the anisotropic 
initial particle momentum distribution using both smeared particles and the NGP method. 
These simulations illustrate the improved accuracy of the smearing method as compared 
to the NGP method. Moreover, thanks to the smaller number of test-particles required by 
the smearing method we are now able to also run simulations on larger lattices. For these 
initial conditions, we observe an instability to transverse magnetic fields while all other field 
amplitudes are essentially constant over the depicted time interval. We also note that the 
evolution of the fields is essentially independent of the lattice spacing, indicating that only 
soft field modes are important for the dynamics of the instability in the Abelian theory. 
Below, we shall find that the situation is different for non-Abelian fields. 

The growing fields deflect the particles if their momenta are not asymptotically large and 
as a consequence they aquire non- vanishing longitudinal momenta Pz (c.f. section [IV]) . We 
remark that PIC methods are not well suited for simulations of the extreme weak-field limit, 
where the particles propagate on nearly straight-line trajectories for a very long time.^ This 
is not due to physical but technical limitations: the weaker the fields, the more accurate 
the particle current J^ needs to be (for example with respect to cancellations of positively 
and negatively charged particles etc.), as it represents a source-term in the field equations. 
High-precision currents, in turn, can only be achieved by employing a large number of 
test particles, driving up the computational time. In practice, the weakest fields for which 
simulations were feasible correspond to initial energy densities of about O.lm^/g'^ for f/(l) 



^ This is better done within the hard-loop approximation. 



10" 



"55 

c 

V 

■a 
0)10== 

V 

c 

LU 



-6^/2 
° 6^/2 
■■■■E2/2 
-^ £2/2 



U(1) 16 Np=18 smearing 




10 



20 



mj 



U(1)64''Np=4 smearing 




10" 



O) 

■T~8 

#10^ 

>. 

c 

V 

■a 
SilO^ 

V 

c 

LU 



— B^/2 
° Bf/2 

■■■■ E^/2 
' E?/2 



U(1) 32^ Np=20 smearing 




30 



10 



20 



mJ 
U(1)16^Np=1000NGP 



30 



i^fi^a^ts 9 msfi fro^a-a^a^-o^ ii^s^fl'*^-*' 4 *xi's <*&***• 



10 20 

mJ 



30 




FIG. 2: Time evolution of the field energy densities for U{1) gauge group and anisotropic particle 
momentum distribution ()2.10p on 16^, 32'^ and 64'^ lattices. Simulation parameters are L = 5 fm, 
Ph = 16 GeV, g"^ Ug = 20 fm-^, moo = 0.1 GeV. 

and ~ 10m^/(7^ for SU{2) gauge group, respectively. The natural regime of applicability 
of PIC methods is when the separation between hard and soft degrees of freedom is not 
asymptotically large. 



C. PIC simulations in non-Abelian gauge theories (CPIC) 

A straightforward extension of the above-mentioned smearing method to the non-Abelian 
case would be to define the current as 



Uhj) = Q^^^{l-Wy), 



j,(z,j + i) = g,^^^iy, 



At 



y ' 



At 



At 



(3.17) 



Jyihj) = Q^^VT^(1 - W,), Jy{i + 1, j) = Q,y-^^W^ , (3.18) 
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where we define the parallel transport of the charge in 2d as ^ 

Q,. = Uiit,j)QU,it,j), Qy ^ Ul{t,j)QUy{t,j) . (3.19) 

One can easily check that this satifies the lattice covariant continuity equation, 

p{t,j) = J2 Ul{i - x).Ut - x)US -x)- -Uhj), (3.20) 

X 

for sites (z, j), (z + 1, j), {i,j + 1): 

pihj) = Jxii,3) + Jyihj), (3.21) 

p(^ + l,j) = Ui{t,j)Ut,j)U,it,j) - Jyit + l,j), (3.22) 

p{t,j + 1) = Ul{t,j)Jy{t,j)Uy{t,j) - UhJ + 1), (3.23) 

Eqs. 03.211) . 03.221) . 03.23P are consistent with the following definitions of the charge densities: 

p{i,j) = Q{l-x){l-y), (3.24) 

p(z,j + l) = Qy{l-x)y , (3.25) 

p(z + l,j) = Q.,x{l-y). (3.26) 

However, since a particle's color charge depends on its path, so does p{i + 1, j + 1) and we 
are not able to calculate it from the charge distribution itself. Rather, we directly employ 
covariant current conservation to determine the increment of charge at site (i + 1, j + 1) 
within the time-step. This way, we can satisfy Gauss's law in the non-Abelian case. The 
forces acting on a particle are determined by analogy to 03.15113.161) . except that the charges 
at neighboring sites is given by the expressions above. 

Finally, we have to check that Tr((5^) is conserved by this smearing method. This is true 
when the lattice spacing a is small, as the total charge of a particle is given by 

Qq = (5(1 - x)(l - y) + Qxx{l -y) + Qy{l - x)y + [apQ^y + (1 - ap)QyJ^xy , (3.27) 

where the Op depend on the path of a particle and Qxy = U^ihj + ^)QyUx{i,j + 1), Qyx = 
U^{i + l,j)QxUy{i + 1, j). If we require that Tr(QQ) be constant, then the cross terms, for 
example Tt{QQx), have to vanish. This is true when a is small, because Tr(Q[A, Q]) = 0: 

Tt{QQx) = Tr(g(g + iga[Ax, Q] + ^(a'))) = Tr{Q^) + ^(a^). (3.28) 

IV. RESULTS 

A. Electrodynamics 

In Fig. [3] we show simulation results for "weak" fields with initial energy density on the 
order of 0.1 m^/g'^. For this case, we find that not only the transverse magnetic but also 



We omit the term [A^, Jz] here, because this term does not appear in 3d and because it does not matter 
for the present discussion. 
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FIG. 3: Time evolution of the field energy densities for U{1) gauge group and anisotropic particle 
momentum distribution on various lattices. Simulation parameters in physical units are L = 5 fm, 
Ph = 8 GeV, ^2 ug = 100 fm-3, ruoo = 0.3 GeV. 

the transverse electric fields grow, albeit at a slower rate.^ Note that here, and throughout 
the manuscript, squared transverse field components denote averages of squares of x- and 
y- components: B^ = {B^ + By)/2, E^ = {E^ + Ey)/2. For the 16^ lattice there is a very 
sudden growth of all field components at t ~ 16m^ which we deem unphysical since we 
do not expect equipartitioning of transverse and longitudinal fields; this is an example for 
the problems which arise in PIC simulations when the fields are very weak and the currents 
are not sufficiently stable^". Indeed, this numerical instability disappears on larger lattices 
with more test-particles. The longitudinal field components are now essentially constant 
and only the transverse components grow. The time where the instability for the transverse 
magnetic fields sets in, as well as their growth rate and saturation point is quite similar for 
the 32^ and 64^ lattices. 

The slope of the exponential growth of B^ is about half that predicted by the hard-loop 
approximation (= y/2moo for the field energy density) which indicates that non-linear effects 



43 



44| 



^ There is little difference to ld-3v simulations, compare for example to Fig. 1 of ref. |42l . 
^'^ Also, on coarse lattices non-linear terms in the Maxwell equations, which arise in the compact lattice 
formulation of U{1), may play a role. This is probably not the case here, however, since qualitatively 
correct results were obtained on the same 16'^ lattice for higher field strength, Fig. [51 
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Final U(1) 64^ 
Final U(1) 32^ 




k^L/Ti 

FIG. 4: Fourier transform of the tranverse magnetic field from Fig. [3l 

are not completely negligible. This is also apparent from the much slower growth of Et as 
compared to B^ and in the Fourier transform of the transverse magnetic field shown in 
Fig. m which illustrates the spectrum of unstable modes. In the linear approximation, the 
spectrum of unstable modes for an extreme particle anisotropy extends to /cmax ~ ''^oo/^, 
where 6 denotes the typical angular width of the particle momentum distribution. For our 
initial condition (12.101) ^ = at t = but grows to ~ 10~^ during the short initial transient 
time before the onset of rapid field growth (see below). However, for the full non-linear 
solution we do not observe such a broad band of instability, modes above k ~ IOtt/L are 
clearly stable. (We note also that the spectrum is not very sensitive to the lattice spacing, 
which indicates that this cutoff is not a lattice artifact.) This might be related to the fact 
that in our simulations 6 increases only little and A;max is not very much smaller than the 
particle scale ph- Therefore, it is likely that back-reaction only allows growth of modes with 



B. Chromodynamics 



We now turn to simulations for the non-Abelian SU{2) gauge group. 3d-3v simulations 
within the hard-loop approximation show that the exponential growth of transverse rnag- 
netic and electric fields saturates at the scale ml^/g due to field self-interactions ^^ [l5|, [l6 . 
While we can not address such weak fields with PIC methods for SU{2), we focus here on 
non-Abelian plasma dynamics beyond the hard-loop approximation for initial field energy 
densities of ~ 10 m'^/g'^ and above. 

Fig. [5] shows results obtained on various lattices. In each case the initial condition was 



^^ This, in fact, is true for moderate anisotropics where the typical longitudinal momentum is less than 
but comparable to the typical transverse momentum. For "extreme" anisotropics such as p.lO|) the field 
strength is expected to grow larger, see discussion below. 
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FIG. 5: Time evolution of the field energy densities for SU{2) gauge group and anisotropic initial 
particle momentum distributions. Simulation parameters are L = 5 fm, Phard = 16 GeV, g^ Ug = 
10/fm3, moo = 0.1 GeV. 



taken to be Gaussian random chromoelectic fields which were "low-pass" filtered such that 
only the lower half of available lattice modes were populated. Note that in Fig. [S] the energy 
densities are plotted on a /mear rather than logarithmic scale. Runs with an isotropic particle 
distribution are shown in each panel as an indication of our numerical accuracy (error bars 
are not indicated for the isotropic runs). The isotropic runs show nearly constant fields over 
the time interval t < 30 m^}^ From these figures we observe that after a period of rapid 
growth both electric and magnetic fields settle to an essentially constant energy density. 

Hence, we confirm that for 3d-3v simulations, even beyond the hard-loop approximation, 
the time evolution of non-Abelian fields stronger than ~ m^/g"^ differs from that in the 
(effectively Abelian) extreme weak-field limit. In particular, a sustained exponential growth 
is absent even during the stage where overall the backreaction on the particles is weak. It 
should be emphasized, in fact, that for very strong anisotropies a linear analysis predicts 
that the exponential field growth (in the weak-field situation) could perhaps continue until 
fi^ ~ m'^/g^/A9'^, with A9'^ = pl/p^ Q. For our initial condition fl230|) . M"^ = at 
t = but grows to (9(10~^) during the initial transient time with constant fields [tm 



< 10 
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In the 64"^ runs the initial amplitudes of the anisotropic and isotropic energy densities are slightly different. 



14 



in Fig. [5]) due to deflection of tfie particles. However, such strong fields are not seen in our 
PIC simulations. This may be related to the above-mentioned back-reaction effects which 
prevent instability of modes near /cmax ~ rrioo/^O. 

The results shown in Fig. \5\ indicate a sensitivity to hard field modes at the ultra-violet 
end of the BriouUin zone, k = 0{a~^), in contrast to the U{1) simulations shown above 



and to earlier ld-3v SU{2) simulations [42, |43|, |4^. The energy density contained in the 



fields at late times increases by a factor of 1.5 when going from a 16^ to 32^ to 64^ lattice 
with the same physical size L. Hence, the dynamics of SU{2) instabilities seen here is not 
dominated entirely by a band of unstable modes in the infrared but clearly involves a cascade 
of energy from those modes to a harder scale A |29|]. However, the simulations shown in 
Fig. [5] indicate that A(t) grows to 0{a~^) during the period of rapid growth of the field 
energy density; otherwise, the final field energy density would not depend on the lattice 
spacing. 

Note that in our PIC simulations the total energy of particles and fields is conserved, 
hence A(t) can not grow arbitrarily large if the energy density of the fields is UV sensitive. 
In three spatial dimensions, this is the case whenever the occupation number of hard field 
modes drops more slowly than ~ l/fc^. The analysis of ref. [29|, for example, suggests a 
~ l/k"^ fall-off near /c ~ A, which indeed corresponds to a quadratically divergent energy 
density as A ^ oo. A different analysis [5^ suggest a ~ 1/k spectrum. 

Although energy conservation will eventually stop the growth of the fields as the lattice 
spacing decreases towards the continuum, it does not solve the following problem. When 
A(t) ~ 1/a, the hard field modes have reached the momentum scale of the particles, ph = 
(9(a~^), and so the clean separation of scales is lost, on which the approach fl2.1|2.2l) is based. 
In fact, since the occupation number (or phase space density) at that scale is of order 1 or 
less by construction^^, it is inappropriate to describe modes at that scale as a classical field. 
Those perturbative modes should be converted dynamically into particles at a lower scale 
A/^p ^ 1/a so that the field energy density and the entire coupled field-particle evolution 
is independent of the artificial lattice spacing. We hope to be able to study this question in 
the future. 

In Fig. O (lower right panel) we have also plotted the tranverse and longitudinal contri- 
butions to the electric and magnetic energy densities as a function of time. From this figure 
we see that at late times the distribution of field energy obtained from our 3d-3v SU{2) 
PIC simulations is approximately isotropic. In contrast, for the ld-3v and 3d-3v U{1) PIC 



simulations and the ld-3v SU{2) CPIC simulations [42|, |43|, |4J| transverse magnetic fields 
dominate, corresponding to an energy-momentum tensor of the form Tg^j^ ^ ^Idd (and 
Afield + ^fiJid ~ 0)- However, we find that with our 3d-3v SU{2) CPIC simulations the energy 
transfered to "hard" field modes is distributed isotropically, regardless of the strong residual 
anisotropy of the particle momenta. The approximate isotropy of the field energy densities 
is also seen in HL simulations of strongly anisotropic plasmas (see below). 

In Fig. Owe show E'^{k) ■ E°-{—k) and B^(k) ■ B"'{—k) (field spectra) as functions of k^ 
at kx = ky = 0. The field configurations at the corresponding time were gauge transformed 
to satisfy Coulomb gauge, diA^ = in continuum notation. The Figure shows that before 
the onset of the rapid field growth from Fig. [5] there is a relatively slow build-up of UV 
modes, in qualitative agreement with the discussion in ref. J29|. Also note the strong growth 



^^ If not, the lattice was chosen too coarse to begin with, and the lattice evolution is far from the continuum 
limit. 
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FIG. 6: Squares of the Fourier transformed (color-) electric and magnetic fields in Coulomb gauge 
(lattice units) at four different times. Upper (lower) panels correspond to 32^ (64^) lattices. Sim- 
ulation parameters are the same as for Fig. Note that tt/L 



IJlr 



of the soft (unstable) field modes at early times. This is the stage where the field evolution 
is independent of the lattice spacing. However, during the period of steep growth from 
Fig. [5] the energy drained from the particles quickly avalanches from IR to UV field modes 
and the distribution over k becomes quite broad. We attribute this avalanche to the UV to 
the self-interaction of the Yang-Mills fields since it is not observed for Abelian f/(l) fields 
(Fig. Hj) even for extreme particle anisotropies. 
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FIG. 7: Broadening of the longitudinal particle momentum distribution for the simulations from 
Fig. [5j The line represents the power-law fit from eq. (j4.2p . 



Fig. [7] depicts the time evolution of the typical longitudinal particle momentum. The 
broadening of the longitudinal distribution during the initial transient trrioo^ 10 is indepen- 
dent of the lattice spacing. During the stage of steep field growth, p^ broadens more rapidly, 
especially for the finer 64^ lattice. Since 



(pD 



dNscatt 2 

dt ^' 



(4.1) 



with qz the typical momentum transfer in ^-direction and dNgcatt/dt the collision rate of 
particles with the field modes, it follows that during this stage either the scattering rate 
and/or the typical momentum transfer increase. This observation (and the fact that the 
effect is much stronger for the 64^ lattice) is consistent with the idea that during this stage 
the UV-cutoff A(t) for the classical field grows to 0{a~^). 

A power-law fit to the late-time evolution of (pi) of the form 



{pD 



{pI + pI) 



ocT 



(4.2) 



yields a ~ 1, corresponding to a random walk of the particles in pz with constant collision 
rate and average momentum transfer. ^^ In an expanding metric, this should then lead to a 
~ r~^/^ drop of the typical pz, as argued by Bodeker J53|] (see also section V in [28|). Hence, 
the anisotropy of the hard modes which develops in the early stage of a heavy-ion collision 
due to the rapid longitudinal expansion should be somewhat smaller than expected within 
the original "bottom-up" scenario [J], which predicted y/{pl) ~ r^^/^. 

For even stronger initial fields, we find that the growth rate of the electric and magnetic 
fields and the dependence on the lattice spacing diminishes. Only a very mild linear growth 



^^ Hence, A appears to be nearly constant in the late stages, which is natural if it has already grown to 
O(a-i). 
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particle momentum. The line represents the power-law fit from eq. ()4.2p . The remaining simulation 
parameters are ph = 16 GeV, g'^ Ug = 10/fm^, L = 5 fm, as before. 



remains over time spans which we can access in our simulations, see Fig. |8l A power-law fit 
of the form (14.21) again shows that the variance of the longitudinal momentum distribution 
grows linearly in time, corresponding to a random walk of the particles in p^ with constant 
scattering rate (off the field modes) and typical momentum transfer. 

We have also performed a simulation with a more moderate initial anisotropy of 
\/(Pz)/(Pt) = 10%, replacing the function 6{pz) in (I2.10p by an exponential ~ 
exp(— 10 \pz\/Ph)- The time evolution of the fields and of their Fourier transform is shown 
in Fig. [91 As expected, the onset of field growth is delayed, by roughly a factor of two, as 
compared to the simulation with extreme particle anisotropy. Also, the intermediate rapid 
growth of the fields is weakened since the rate of energy transfer from particles to fields 
is now lower. Nevertheless, even for this run with more moderate anisotropy, we observe 
that the Fourier spectrum flattens somewhat at late times (when the growth saturates); it 
actually resembles the final spectrum from Fig. O although the angular width of the particle 
momentum distribution is now much larger. 



Comparison with hard-loop results 

In Fig. [TU]we present the dependence of the energy densities obtained by solving the soft- 
field equations of motion in the hard-loop approximation 16| for the case of ^/{pI)/{Pt) ~ 
0.1. For this simulation the initial condition has no fields, but instead has finite fluctuations 
of the auxilliary W fields which induce field fluctuations with typical initial energy density 
S ~ 0.5 m^/g'^ at early times (for details of the initial condition and numerical algorithm 
used see Ref. Il6|). Note that if one makes the same run with an isotropic hard particle 
distribution the fields remain isotropic and show no signs of growth. 

As can be seen from Fig.dOlthere is an initial period where no field growth occurs, similar 
to Fig. [9l At TUoot ~ 10 all field components begin to grow isotropically in the hard- loop 
simulation. They then grow approximately linearly and isotropically for the duration of 
the run. At late times one can see a slight curvature upwards suggesting a faster than 
linear growth. This curvature is related to the finite discretizations used in both space-time 
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FIG. 9: Field evolution for moderate initial anisotropy, \ (pi) / ipj) = 10%. The remaining simula- 
tion parameters are the same as in Fig. [3 Top: field energy density as a function of time. Bottom: 
Fourier transformed color-fields (in Coulomb gauge) at different times. 



and velocity space. Improving these discretizations or damping high angular momentum 
modes [l5|] extends the region in time over which the linear behavior is observed. The 
inset of Fig. [10] shows the dependence of the field energy densities at late times when the 
magnitudes are comparable to those shown in Fig. [91^^ Also, similar to the CPIC results 
we see that the energy densities obtained in the hard-loop simulation are approximately 
isotropic. 



^^ Note that to match to the total energy density take into account the factor of three difference between 
the total and the contributions from individual components. 
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FIG. 10: Field evolution for moderate initial anisotropy as obtained from numerical solution of 
hard- loop effective theory 16]. 



V. DISCUSSION 



We have presented numerical solutions of coupled particle-field dynamics in Abelian and 
non-Abelian gauge theories beyond the HL approximation. Such plasmas exhibit collective 
instabilities when the distribution of particles in momentum space is anisotropic (we restrict 
ourselves to an oblate distribution here). However, in the regime where field self-interactions 



can not be neglected, corresponding to field strengths of order 



~ m: 



Jg and above, the evolu- 



tion in chromodynamics is rather different from that in Abelian plasmas. Most remarkably, 
for the initial conditions employed here, we do not observe a period of sustained exponential 
growth of the soft fields. This is because the field modes can interact directly and so the 
energy drained by the instability from the hard particles does not remain in unstable modes 
for a sufficiently long time. 

The transition from exponential to linear growth at the scale rri^/g has been seen previ- 



ously in 3d-3v simulations within the HL approximation [15|, [iGj. It appears that if the par- 
ticles exhibit only weak anisotropy, that the HL approach could be continued until very late 
times since hard field modes are populated very slowly. However, for moderate anisotropics 
the instability is also affected more easily by collisions among the particles [54 and by the 
rapid longitudinal expansion encountered in relativistic heavy- ion collisions [l7|. 

For very anisotropic particle momentum distributions we find here that the evolution 
differs significantly for field strengths of order rn^/g and above, even if overall the back- 
reaction on the hard particle modes is small. After some initial transient time, we observe 
a short period of rapid growth of soft electric and magnetic fields followed by an avalanche 
of energy into the UV. As long as the energy of the classical Yang-Mills fields is still much 
less than that of the particles, this avalanche proceeds to the smallest length scales, namely 
the lattice spacing a, before reaching a steady-state distribution. Since 1/a will eventually 
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reach ph as one takes the continuum hmit, the loss of a clear separation of scales between 
soft and hard degrees of freedom could perhaps be viewed as (partial) "classical decay" of 
the field into particles. To understand what happens then requires us to go beyond the 
Wong- Yang-Mills equations, which rely on such a separation of scales. We also point out 
that the energy-momentum tensor of the fields is nearly isotropic, while Weibel instabilities 
seen in U{1) and ld-3v SU{2) simulations lead to dominance of transverse magnetic fields. 
If the initial field energy density is rather high, of order (10^ — 10"^) m'^/g'^, we do not 
observe an instability with rapid growth of the fields. In this case, the mean-square width 
of the longitudinal particle distribution grows linearly in time. This corresponds to a ran- 
dom walk with constant scattering rate (from the field fluctuations) and typical momentum 
transfer. For even stronger initial fields with energy density comparable to that of the par- 
ticles, ~ PhU n, w e of course expect a very rapid isotropization of the particle pressure, see 



e.g. refs. [42, |43|, |44 



Our simulations indicate that it could be important to implement a dynamical field to 
particle conversion at a scale Aj^p < ph to avoid that field and particle degrees of freedom 
overlap. This might lead to a continuous "cycle" of energy transfer from the anisotropic 
particles to soft field modes which then quickly re-emerges at that scale in isotropic form. 
The cycle should continue until all particle momenta are isotropic. Theoretically, it also 
remains to be understood why the energy avalanche from soft to hard field modes occurs 
so rapidly. An independent analytical discussion of energy transfer to small length scales 
in QCD has been recently given in ref. 52]. In the future, we intend to perform more 



quantitative studies and verify the matching to the "bottom-up" solution for particle-field 



thermalization [55|. To do so, however, requires us to also allow for longitudinal expansion 



of the lattice and to include hard elastic and inelastic scattering of particles. 
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APPENDIX A: MOMENTUM UPDATE 

We consider the discretization of the equation 

'^ = gQ\E'' + vxB''). (Al) 

(JjV 

The difference form of the equation is 

p{t + At/2) - p{t - At/2) ^J^a,,, , p(t + At/2) + p(t - At/2) ^ B%t)\ 
= gQ^Eit) + x-^J (A2) 
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Eq. (lA2p can be solved by the Buneman-Boris method 45|, |46j as follows: 



p{t) = p{t-At/2) + ^Eit), (A3) 

e{t) = \pit)\, (A4) 

P'it) =pW + ^pWx^' (A5) 

2 At ,, , B(t) ,,^, 

P^W = '''*' + l + (B/eWAi/2)-^ TP'" ^ 7(77 • '^" 

p(i + Ai/2) = P2(i) + ^-B(«). (A7) 

where E = gQ"-E°- and -B = gQ^-B"-. This scheme is time reversible and the overall 
momentum integration is accurate to second-order in the time step. 
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